/*************************************************************
* m3d_noise.c
* Noise support class implementation. 
*
* Source book: "Texturing and Modeling", David S. Ebert, 
*               F. Kenton Musgrave, Darwyn Peachey, Steven Worley, 
*               Ken Perlin
* Prefix   : m3dNoise
*************************************************************/
#include <stdafx.h>

#include <stdlib.h>

//#include <ap_comm.h>
#include <m3d\m3d.h>

//
// Global variables
//
m3dNOISE      *m3dNoise;

//
// Local Defines
//
#define _M3D_N_INTERP_LAGR_KP  4

#define _M3D_INDEX_CLAMP(i)   ((i)&(size-1))
#define _M3D_INDEX_HASH(i)    (tabHash[_M3D_INDEX_CLAMP(i)])
#define _M3D_INDEX(ix,iy,iz)  _M3D_INDEX_HASH((ix) + _M3D_INDEX_HASH((iy) + _M3D_INDEX_HASH((iz))))

//
// Static variables
//
static m3dSPL_LAGRANGE  *_m3dSplLagr;

//
// Static functions
//
static float _m3dSplLagrGetInterpValue (float *y, float arg);
static BOOL  _noiseFillHash            (int *tabHash, int size);

/********************************************************************
* m3dInitNoise ()
* Initialize noise subsystem
* IN : none
* OUT: none
* RET: FALSE if error occured
********************************************************************/
BOOL m3dInitNoise (void)
{
   int i;

   if ((_m3dSplLagr = apNEW m3dSPL_LAGRANGE(_M3D_N_INTERP_LAGR_KP)) == NULL) {
      return FALSE;
   }
   for (i=0; i<_M3D_N_INTERP_LAGR_KP; i++) {
      _m3dSplLagr->SetKpArg(i, (float)i);
   }

   // create system noise
   if ((m3dNoise = apNEW m3dNOISE_GRAD(512)) == NULL) {
      return FALSE;
   }
   m3dNoise->Fill(0);

   return TRUE;
}

/********************************************************************
* m3dTermNoise ()
* Terminate noise subsystem
* IN : none
* OUT: none
* RET: none
********************************************************************/
void m3dTermNoise (void)
{
   delete _m3dSplLagr;
   delete m3dNoise;

   return;
}
/*********************************************************************
* m3dNoiseCatRom2()
* IN : x - argument
* OUT: none
* RET: value of the Catmull-Rom spline
* REFERENCE: [1], p.78-79
*********************************************************************/
__INLINE float m3dNoiseCatRom2(float x)
{
   float catrom2;

   if (x > 4) {
      return 0.f;
   }

   if (x < 1.f) {
      catrom2 = 0.5f * (2.f + x*x*(-5.f + 3*x));
   } else {
      catrom2 = 0.5f * (4.f + x*(-8.f + x*(5.f - x)));
   }
   return catrom2;
}

/*************************************************************
* m3dNOISE Class Members
*************************************************************/
/*************************************************************
* m3dNOISE::Fill ()
* Init noise object
* IN:  type - type of noise (NOISE_TYPE_...)
*      size - size of noise table (quality)
* OUT: none
* RET: pNoise or FALSE if error occuerd
*************************************************************/
BOOL m3dNOISE::Fill (int seed)
{
   if ((tabHash = (int*)realloc(tabHash, size*sizeof(int))) == NULL) {
      return FALSE;
   }
   if ((table = realloc(table, size * cellSize)) == NULL) {
      return FALSE;
   }
   if (!_noiseFillHash(tabHash, size)) {
      return FALSE;
   }
   return TRUE;
}
/*************************************************************
* m3dNOISE::~m3dNOISE ()
*
* Noise object destructor
* IN:  none
* OUT: none
* RET: none
*************************************************************/
m3dNOISE::~m3dNOISE (void)
{
   free(table);
   free(tabHash);
   return;
}
/*************************************************************
* m3dNOISE_VALUE_PROTO::Fill ()
*     Fill value noise table
* IN : seed   - randomization seed
* OUT: none
* RET: FALSE if an error occured
*************************************************************/
BOOL m3dNOISE_VALUE_PROTO::Fill (int seed)
{
   float *pTable;
   int    i;

   if (!m3dNOISE::Fill(seed)) {
      return FALSE;
   }

   pTable = (float*)table;
   for (i = 0; i < size; i++) {
      *pTable = 2 * m3dRandom() - 1.f;
      pTable++;
   }
   return TRUE;
}

/*************************************************************
* m3dNOISE_GRAD::Fill ()
*     Fill gradient noise table
* IN : seed   - randomization seed
* OUT: none
* RET: FALSE if an error occured
* REFERENCE: [1], pp.72-76
*************************************************************/
BOOL m3dNOISE_GRAD::Fill (int seed)
{
   float  z, r, theta;
   float *pTable;
   int    i;

   if (!m3dNOISE::Fill(seed)) {
      return FALSE;
   }

   pTable = (float*)table;
   for (i = 0; i < size; i++) {
      // produce uniformely distributed vector on the unit sphere
      // get it in the polar coords system:
      // z     - sine of latitude angle of vector,
      // r     - cosine of latitude angle,
      // theta - longitude angle of vector
      z = 1.f - 2.f * m3dRandom();
      r = (float)m3dSqrt(1 - m3dPow2(z));
      theta = 360.f * m3dRandom();
      // point of the (x, y) circle
      pTable[0] = r * m3dCos(theta);
      pTable[1] = r * m3dSin(theta);
      pTable[2] = z;
      pTable += 3;
   }
   return TRUE;
}
/*********************************************************************
* m3dNOISE_GRAD::GetValue ()
*     Get gradient noise value
* IN:  position
* OUT: none
* RET: value
*********************************************************************/
float m3dNOISE_GRAD::GetValue (m3dV *pPos)
{
   m3dVINT posI;
   m3dV    posF, w, vDir;
   m3dV   *pValue3D;
   float   v_xyz[2][2][2];
   float   v_xy[2][2];
   float   v_x[2];
   int     j[3], i;
   float   rangeMin, rangeMax;
   float   value;

   for (i=0; i<3; i++) {
      m3dVIntComponent(&posI,i) = m3dIntFloor(m3dVectComponent(pPos,i));
      m3dVectComponent(&posF,i) = m3dVectComponent(pPos,i) - (float)m3dVIntComponent(&posI,i);
      m3dVectComponent(&w,i)    = m3dSmoothStep(m3dVectComponent(&posF,i));
   }

   // interpolate between collinear axis in 3d
   for (j[0]=0; j[0]<2; j[0]++) {
      m3dVectComponent(&vDir,0) = m3dVectComponent(&posF,0) - (float)j[0];
      for (j[1]=0; j[1]<2; j[1]++) {
         m3dVectComponent(&vDir,1) = m3dVectComponent(&posF,1) - (float)j[1];
         for (j[2]=0; j[2]<2; j[2]++) {
            m3dVectComponent(&vDir,2) = m3dVectComponent(&posF,2) - (float)j[2];
            pValue3D = (m3dV*)(*this)[_M3D_INDEX(posI.x+j[0], posI.y+j[1], posI.z+j[2])];
            v_xyz[j[0]][j[1]][j[2]] = m3dDotProduct(pValue3D, &vDir);
         }
         v_xy[j[0]][j[1]] = _m3dLerp(v_xyz[j[0]][j[1]][0], v_xyz[j[0]][j[1]][1], w.z);
      }
      v_x[j[0]] = _m3dLerp(v_xy[j[0]][0], v_xy[j[0]][1], w.y);
   }
   value = _m3dLerp(v_x[0], v_x[1], w.x);

   // make gradient noise more uniformly distributed across the spectrum
   // most of the noise power is in range [-0.5, 0.5], but histogramm
   // could be analized for better spectrum distribution
   rangeMin = -0.5f;
   rangeMax =  0.5f;
   // get rid of the other part of the spectrum
   value  = m3dClamp(-1.f, 1.f, rangeMin, rangeMax, value);

   return value;
}

/*********************************************************************
* m3dNOISE_VALUE::GetValue ()
*     Get value noise value
* IN:  position
* OUT: none
* RET: value
*********************************************************************/
float m3dNOISE_VALUE::GetValue (m3dV *pPos)
{
   m3dVINT  posI;
   m3dV     posF;
   int      j[3], i;
   float    v_xyz[_M3D_N_INTERP_LAGR_KP][_M3D_N_INTERP_LAGR_KP][_M3D_N_INTERP_LAGR_KP];
   float    v_xy[_M3D_N_INTERP_LAGR_KP][_M3D_N_INTERP_LAGR_KP];
   float    v_x[_M3D_N_INTERP_LAGR_KP];
   float    tmp[_M3D_N_INTERP_LAGR_KP];
   float    shift = (float)(_M3D_N_INTERP_LAGR_KP/2);
   float    value;

   for (i=0; i<3; i++) {
      m3dVIntComponent(&posI,i) = m3dIntFloor(m3dVectComponent(pPos,i));
      m3dVectComponent(&posF,i) = m3dVectComponent(pPos,i) - (float)m3dVIntComponent(&posI,i);
   }

   for (j[0]=0; j[0]<_M3D_N_INTERP_LAGR_KP; j[0]++) {
      for (j[1]=0; j[1]<_M3D_N_INTERP_LAGR_KP; j[1]++) {
         for (j[2]=0; j[2]<_M3D_N_INTERP_LAGR_KP; j[2]++) {
            v_xyz[j[0]][j[1]][j[2]] = *((float*)(*this)[_M3D_INDEX(posI.x + j[0]-_M3D_N_INTERP_LAGR_KP/2, 
                                                                     posI.y + j[1]-_M3D_N_INTERP_LAGR_KP/2, 
                                                                     posI.z + j[2]-_M3D_N_INTERP_LAGR_KP/2)]);
         }
         for (i=0; i<_M3D_N_INTERP_LAGR_KP; i++) {
            tmp[i] = v_xyz[j[0]][j[1]][i];
         }
         v_xy[j[0]][j[1]] = _m3dSplLagrGetInterpValue(tmp, posF.z+shift);
      }
      for (i=0; i<_M3D_N_INTERP_LAGR_KP; i++) {
         tmp[i] = v_xy[j[0]][i];
      }
      v_x[j[0]] = _m3dSplLagrGetInterpValue(tmp, posF.y+shift);
   }
   value = _m3dSplLagrGetInterpValue(v_x, posF.x+shift);
   return value;
}

/*********************************************************************
* m3dNOISE_VALUE_CONV::GetValue ()
*     Get value-convolution noise value
* IN:  position
* OUT: none
* RET: value
*********************************************************************/
float m3dNOISE_VALUE_CONV::GetValue (m3dV *pPos)
{
   int      j[3], i;
   m3dVINT  posI;
   m3dV     posF, d;
   float    value, sum = 0.f;

   // detect point's position of base cell 
   // and relative position in the associated cell.
   for (i=0; i<3; i++) {
      m3dVIntComponent(&posI,i)  = m3dIntFloor(m3dVectComponent(pPos,i));
      m3dVectComponent(&posF,i) = m3dVectComponent(pPos,i) - m3dVIntComponent(&posI,i);
   }

   for (j[2] = -2; j[2] <= 2; j[2]++)  {  
      d.z = j[2] - posF.z;
      d.z = m3dPow2(d.z);
      for (j[1] = -2; j[1] <= 2; j[1]++)  {
         d.y = j[1] - posF.y;
         d.y = m3dPow2(d.y);
         for (j[0] = -2; j[0] <= 2; j[0]++)  {
            d.x = j[0] - posF.x;
            d.x = m3dPow2(d.x);
            value = *((float*)(*this)[_M3D_INDEX(posI.x + j[0], posI.y + j[1], posI.z + j[2])]);
            sum += value * m3dNoiseCatRom2(d.x + d.y + d.z);
         }
      }
   }
   return sum;
}

/*************************************************************
* m3dNOISE_SPARSE_CONV::Fill ()
*     Fill sparse convolution noise table
* IN : seed   - randomization seed
* OUT: none
* RET: FALSE if an error occured
*************************************************************/
BOOL m3dNOISE_SPARSE_CONV::Fill (int seed)
{
   float *pTable;
   int    i;

   if (!m3dNOISE::Fill(seed)) {
      return FALSE;
   }
   
   pTable = (float*)table;
   for (i = 0; i < size; i++) {
      // produce uniformely distributed vector on the unit sphere
      // get it in the polar coords system:
      // z     - sine of latitude angle of vector,
      // r     - cosine of latitude angle,
      // theta - longitude angle of vector
      pTable[0] = m3dRandom();
      pTable[1] = m3dRandom();
      pTable[2] = m3dRandom();
      pTable[3] = 1.f - 2*m3dRandom();
      pTable += 4;
   }
   return TRUE;
}
/*********************************************************************
* m3dNOISE_SPARSE_CONV::GetValue ()
*     Get sparse-convolution noise value
* IN:  position
* OUT: none
* RET: value
*********************************************************************/
float m3dNOISE_SPARSE_CONV::GetValue (m3dV *pPos)
{
   int     i, h, n;
   int     j[3];
   m3dVINT posI;
   m3dV    posF, d;
   float   *fp, sum = 0.f;

   // detect point's position of base cell 
   // and relative position in the associated cell.
   for (i=0; i<3; i++) {
      m3dVIntComponent(&posI,i) = m3dIntFloor(m3dVectComponent(pPos,i));
      m3dVectComponent(&posF,i) = m3dVectComponent(pPos,i) - (float)m3dVIntComponent(&posI,i);
   }

   // for all cells at the radius (theoretically - 2, 
   // but since we aren't moving at the 0X & 0Z directions - so, not so 
   // accuracy may be performed)
   for (j[0] = -2; j[0] <= 2; j[0]++)  {
      for (j[1] = -2; j[1] <= 2; j[1]++)  {
         for (j[2] = -2; j[2] <= 2; j[2]++)  {
         // at the 0Y direction we are moving, so apply full range filtering
         // for (j = -1; j <= 1; j++)  {
            h = _M3D_INDEX(posI.x + j[0], posI.y + j[1], posI.z + j[2]);
            // two impulses uniformly distributed at each cell...
            for (n = 1; n > 0; n--, h = (h+1)%size)  {
               fp = (float*)(*this)[h];
               for (i=0; i<3; i++) {
                  m3dVectComponent(&d,i) = m3dVectComponent(&posF,i) - (j[i] + *fp);
                  fp++;
               }
               sum += m3dNoiseCatRom2(m3dLengthVector_2(&d)) * (*fp);
            }
         }
      }
   }
   return sum;
}

/******************************************************************
* Static functions
******************************************************************/
/******************************************************************
* IN  : y - array of the length _M3D_N_INTERP_LAGR_KP
******************************************************************/
static float _m3dSplLagrGetInterpValue(float *y, float arg)
{
   float  result;
   int    i;

   for (i = 0; i < _M3D_N_INTERP_LAGR_KP; i ++) {
      _m3dSplLagr->SetKpValue(i, &y[i]);
   }
   _m3dSplLagr->GetValue(arg, &result);
   return result;
}

#ifdef _AP_XENON
#pragma optimize("", off) // temporary hack for XDK:7978.3
#endif

/*************************************************************
* _noiseFillHash()
* Generate permutation table for the noise
* IN : tabHash - pre-allocated array of 'size' integers
*      size    - table size
* OUT: fill tabHash with random permutation of the set {0,1,2,...,size-1}
* RET: none
*************************************************************/
static BOOL _noiseFillHash(int *tabHash, int size)
{
   int  ind, indHash, indWant;
   int *tabInd;

   // create/init Erastofen's array
   ASSERT(size > 0);
   tabInd = (int*)calloc(size, sizeof(int));
   if (!tabInd) {
      ASSERT(FALSE);
      return FALSE;
   }
   for (ind = 0; ind < size; ind++)  {
      tabInd[ind] = ind;
   }

   for (indHash = 0; indHash < size; indHash++)  {
      // find unused item with given RND position
      ind = -1;
      indWant = (int)(m3dRandom() * (size - indHash));
      STRONG_ASSERT(indWant < (size - indHash));
      do {
         // skip unused elements
         do   {
            ind++;
            STRONG_ASSERT(ind < size);
         } while (tabInd[ind] == -1);

         indWant--;
      } while (indWant > 0);

      // store it!
      tabHash[indHash] = tabInd[ind];
      tabInd[ind]  = -1;
   }
   free(tabInd);
   return TRUE;
}

#ifdef _AP_XENON
#pragma optimize("", on)
#endif

//
// end of file m3d_noise.C
//
